In [3]:
# Some boilerplate code
from __future__ import division
import numpy as np
import pandas
import matplotlib.pyplot as plt
%matplotlib inline
%cd data/
cytosol,advection = None,None
plt.rcParams['figure.figsize'] = 15,7
In [10]:
def get_P_param(**kargs):
param = {}
#Discretisation
param["L"] = 134.6 #µm
param["T"] = 700 #s
# Kinetics
param["kon"] = 4.74e-2 # μm s−1
param["koff"] = 7.3e-3 # s-1
param["psi"] = 0.174 #µm-1
# Concentration
param["rho"] = 1 #µm--3
# Diffusion
param["D"] = 0.15 #µm^2.s-1
# Advection parameter
param["xi"] = 1 #adimensional
#interaction
param["alpha"] = 1
param["beta"] = 2
param["kap"] = 0.19
param["kpa"] = 0.01
for k,v in kargs.items():
param[k] = v
return param
def get_A_param(**kargs):
param = get_P_param()
param["kon"] = 8.58e-3 # µm s-1
param["koff"] = 5.4e-3 # s-1
param["xi"] = 1 # adimensional
param["rho"] = 1.56 #µm--3
param["D"] = 0.28
for k,v in kargs.items():
param[k] = v
return param
# If you want to get a new set of parameters just call:
theta = get_P_param()
# You can specifie the value of some of the parameters:
theta = get_A_param(a=4,D=1000)
In [16]:
def artificial_flow(Nt,Nx,s=1e-2):
flow_artificial = np.zeros((Nt, Nx)) - s
flow_artificial[:,int(Nx/2):] = -flow_artificial[:,int(Nx/2):]
flow_artificial[:,:int(0.2*Nx)] = 0
flow_artificial[:,-int(0.2*Nx):] = 0
flow_artificial[:,int(Nx/2)-0.05*int(Nx/2):int(Nx/2)+0.05*int(Nx/2)] = 0
return flow_artificial
## Curated public version of the notebook: upublished data are not yet avaiable.
#flow = pandas.read_csv("flow-myosin-um-per-min_3040linear-itp",sep=None,header=None).values / 60
#ppar = pandas.read_table("pPAR.txt",header=None,sep=None).values
#apar = pandas.read_table("aPAR.txt",header=None,sep=None).values
## end of curated part
flow = artificial_flow(3000,150)
In [12]:
def diffusion(p,param,dt,l,t):
return (param["D"] * l**-2 *
(np.concatenate((p[1:],[p[0]]))
- 2*p[:]
+ (np.concatenate(([p[-1]],p[:-1])))))
def cytosol(p,param,dt,l,t):
return param["cyto"][t]*param["kon"]*param["psi"] - p*param["koff"]
def mechano(p,a,param,param_a,dt,l,t):
dp = np.zeros(p.shape)
da = np.zeros(a.shape)
param["F"][t+1] = param["F"][t] + param["Fp"]*p + param_a["Fa"]*a
param_a["F"][t+1] = param["F"][t]
return dp,da,param,param_a
def advection(p,param,dt,l,t):
return ( np.concatenate(([param["F"][t,-1]],param["F"][t,:-1])) * l**-1 * param["xi"] * np.concatenate(([p[-1]],p[:-1]))
- np.concatenate((param["F"][t,1:],[param["F"][t,0]])) * l**-1 * param["xi"] * np.concatenate((p[1:],[p[0]]))
)
def interaction(p,a,param,param_a,dt,l,t):
dp = - param["kpa"] * p * a**param["alpha"]
da = - param_a["kap"] * p**param_a["beta"] * a
param["cyto"][t+1] = param["cyto"][t] - dt*np.sum(dp)
param_a["cyto"][t+1] = param_a["cyto"][t] - dt*np.sum(da)
return dp,da,param,param_a
In [13]:
def play2(p,a,param,param_a,components,components2=[]):
""" The simplest form of the model: Diffusion only
Args:
p (np.array): A timexSpace array, the first line give
the initial condtion, its shape gives the number of
discretization volumes.
a (np.array): A timexSpace array, the first line give
the initial condtion, its shape gives the number of
discretization volumes.
param (dict): Parameters for p
param_a (dict): Parameters for a
components (tuple): the list of functions to use in the computation of dP/dA.
components_2 (tuple): the list of functions to use in the computation of dP&dA.
Returns:
(np.array): with a timestep by line.
"""
dt = param["T"]/float(p.shape[0]-1) #Time discretisation.
l = param["L"]/float(p.shape[1]-1) #Space discretisation.
param["dt"] = dt
param["l"] = l
print "Simulation of 2 species with {},{}, T={} (dt={}), L={} (dx={})".format([x.__name__ for x in components],
[x.__name__ for x in components2],
param["T"],dt,param["L"],l)
#Cytosol
param["cyto"],param_a["cyto"] = np.zeros(p.shape[0]), np.zeros(a.shape[0])
param["cyto"][0] = param["rho"] - param["psi"]/p.shape[1] * np.sum(p[0,:])
param_a["cyto"][0] = param_a["rho"] - param_a["psi"]/a.shape[1] * np.sum(a[0,:])
# Checking for advection stability
if advection in components:
dmax = max(param["F"].max() * param["dt"],param_a["F"].max() * param["dt"])
if dmax >= 0.1*l:
print ("Warning ! Maximal flow movment in one step is {}µm "
"whereas the space discretization length is {}µm.\n "
"You could try to reduce dt to improve that.").format(dmax,l)
# Time loop
for t in range(p.shape[0]-1):
dpoverdt = []
daoverdt = []
for fct in components:
dpoverdt.append(fct(p[t,:],param,dt,l,t))
daoverdt.append(fct(a[t,:],param_a,dt,l,t))
if cytosol in components:
param["cyto"][t+1] = param["cyto"][t] - dt*np.sum(cytosol(p[t,:],param,dt,l,t))
param_a["cyto"][t+1] = param_a["cyto"][t] - dt*np.sum(cytosol(a[t,:],param_a,dt,l,t))
for fct in components2:
x,y,param,param_a = fct(p[t,:],a[t,:],param,param_a,dt,l,t)
dpoverdt.append(x)
daoverdt.append(y)
p[t+1,:] = (p[t,:] + dt*(np.sum(dpoverdt,0)))
a[t+1,:] = (a[t,:] + dt*(np.sum(daoverdt,0)))
return p,a,param,param_a
In [14]:
def display2(p,a,param,param_a,cyto=False,flow=False):
"""Display function
Args:
p (np.array): A timeXspace array.
param (dict): the parameters
cyto (bool): display the membrane/cytosol concentration"""
vmin = min(p.min(),a.min())
vmax = max(p.max(),a.max())
plt.subplot(1,2,1)
plt.imshow(p,interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0),
vmin=vmin,vmax=vmax)
plt.title("P")
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(a,interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0),
vmin=vmin,vmax=vmax)
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.title("A")
plt.colorbar()
if cyto:
plt.show()
plt.subplot(2,1,1)
plt.plot(param["cyto"],label="P Cytosol")
plt.plot(p.sum(1),label="P Membrane")
plt.plot(p.sum(1)+param["cyto"],label="P Overall")
#print max(p.sum(1)+param["cyto"]),min(p.sum(1)+param["cyto"])
plt.xlabel("time (sec)")
plt.ylabel("Concentration")
plt.legend()
plt.subplot(2,1,2)
plt.plot(param_a["cyto"],label="A Cytosol")
plt.plot(a.sum(1),label="A Membrane")
plt.plot(a.sum(1)+param_a["cyto"],label="A Overall")
plt.xlabel("time (sec)")
plt.ylabel("Concentration")
plt.legend()
if flow:
plt.show()
plt.subplot(1,2,2)
plt.imshow(param["F"],interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0),
vmin=vmin,vmax=vmax)
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.title("Flowspeed um.sec-1")
In [18]:
P = np.zeros(flow.shape)
A = np.zeros(flow.shape)
PARAMETERS_P = get_P_param(F=flow,)
PARAMETERS_A = get_A_param(F=flow,)
results = play2(P,A, # Grid
PARAMETERS_P,PARAMETERS_A, # Parameter
[diffusion,cytosol,advection],# First order phenomenon
[]) # Second order phenomenon
display2(*results,cyto=True)
plt.show()
## Curated public version of the notebook: upublished data are not yet avaiable.
#print "Experimental results:"
#display2(ppar,apar,get_P_param(),get_A_param())
## end of curated part